NASA 

Technical 

Paper 

2688 

May 1987 


Quantitative Analysis 
of the Reconstruction 
Performance of Interpolant 


Donald L. Lansing 
and Stephen K. Park 


NASA 


NASA 

Technical 

Paper 

2688 

1987 


Quantitative Analysis 
of the Reconstruction 
Performance of Interpolant 


Donald L. Lansing 

Langley Research Center 
Hampton, Virginia 

Stephen K. Park 

The College of William and Mary 
Williamsburg, Virginia 


NASA 

National Aeronautics 
and Space Administration 


Scientific and Technical 
Information Office 


Contents 


Introduction 1 

Symbols 1 

Statement of the Problem 1 

Analysis 2 

Linear, Shift-Invariant Interpolants 3 

Mean-Square-Error Criterion 3 

Interpolation Methods Investigated 4 

Local Splines ‘ 4 

Parametric cubic Hermite splines 4 

Parametric quintic Hermite splines 5 

Global Splines 5 

Exponential splines 5 

Nu splines 5 

Global interpolation functions 5 

Reconstruction Filters 6 

Parametric cubic convolution (PCC) 6 

Keys’ cubic 6 

BAWA cubic 6 

Quadratic Shape-Preserving Splines 6 

Comparison of Reconstruction Properties of Interpolants 7 

Interpolation Function 7 

Reconstruction Filter 8 

The Function e 2 (v) 8 

Optimized Quintic Hermite Splines 8 

Optimized Parametric Cubic Interpolants 9 

Concluding Remarks 10 

Appendix A — Equations for Hermite Basis Splines 12 

Parametric Cubic Hermite Splines 12 

Parametric Quintic Hermite Splines 12 

Appendix B — Interpolation Functions 13 

Appendix C — Some Analytic Expressions for f(i /) and e 2 (v) 15 

References 16 

Table I 17 

Figures 18 


PRECEDING PAGE BLANK NOT nJ0& 

* 


Introduction 


Interpolation is a technique of fundamental im- 
portance in many disciplines including signal and im- 
age processing, computer graphics, computer-aided 
design (CAD), and numerical analysis. In all these 
disciplines, one is frequently faced with some ver- 
sion of the following problem: given a discrete set 
of points which correspond to values of a function 
y(x) sampled on a grid, construct an easily com- 
puted function g(x) which equals y(x) on the sam- 
pling grid and approximates y(x) elsewhere. In the 
signal- and image-processing literature, g(x) is said 
to reconstruct y(x) from its samples; in the computer 
graphics, CAD, and numerical analysis literature the 
term interpolate is more common. In either case, the 
concept is the same. 

There are numerous interpolation and reconstruc- 
tion algorithms from which to choose. Most of these 
algorithms are based upon interpolants which are ei- 
ther polynomial or exponential splines. If the number 
of samples is small, global algorithms are generally 
preferred. However, if the number of samples is large, 
the algorithms are virtually always local. Generally, 
global splines are used in CAD, graphics, and nu- 
merical analysis applications, whereas local splines 
are employed for signal- and image-processing prob- 
lems. In either case, choosing an interpolant for a 
specific application involves a combination of subjec- 
tive factors including ease of use and the analyst’s 
prior experience and personal judgment as to what 
constitutes a “good” technique. 

This paper is concerned with the application of 
a quantitative measure of the performance of vari- 
ous interpolation and reconstruction algorithms. The 
analysis provides a quantitative criterion based upon 
the mean square error for assessing how well g{x) 
approximates y(x). It is applicable to the situation 
commonly occurring in signal- and image-processing 
in which the number of grid points is very large and 
the grid is equally spaced. The interpolants under 
consideration are characterized by the properties of 
linearity and shift invariance; that is, g(x) is a lin- 
ear combination of the sample values and a left or 
right shift of the sample grid produces only a corre- 
sponding shift in g(x). The criterion, reformulated 
in the frequency domain, is applicable to data sets 
with arbitrary spectral content. However, this pa- 
per is primarily concerned with the interpolation of 
sufficiently sampled, band-limited data. The results 
provide a unified approach which applies equally well 
to evaluating interpolants for CAD, graphics, numer- 
ical analysis, and image-processing applications. 


Symbols 


e 2 (u) 

function defined by the series 
in equations (6), expresses an 
interpolant ’s contribution to the 
mean square error 

0(x) 

interpolating function 

Hi(s) 

quintic Hermite basis functions, 
* = 1, 2, 3, ..., 6 

h{(s) 

cubic Hermite basis functions, 
i = 1, 2, 3, 4 

i, n 

knots, sample grid locations 

K 

integer defining width of interpola- 
tion function 

r{x) 

interpolation function 

r(v) 

reconstruction filter, Fourier trans- 
form of r(x) 

s 

= x — 2 , independent variable in 
Hermite function 

x 

independent variable of sampled 
function 

y(n) 

data samples of y(x) at the knots 

y{x) 

function to be sampled and 
interpolated 

yW 

Fourier transform of y(x) 

\yH\ 2 

data energy spectrum 

a, (3, t, t' 

parameters in Hermite interpolation 
and parametric cubic convolution 

€ 2 

average mean square error 

V 

frequency parameter, cycles per 
sample 

Vc 

cut-off frequency 

G 

shape parameter in model of data 
energy spectrum 

Abbreviations: 

BAWA 

piecewise cubic interpolant 
developed in reference 9 

CAD 

computer-aided design 

PCC 

parametric cubic convolution 

OBE 

unfiltered out-of-band energy 


Statement of the Problem 

This paper is concerned with the following ver- 
sion of the standard interpolation problem. Sup- 
pose that the real- valued function y(x), defined 


for all real x, has been sampled at the integers 
—2, — 1, 0, 1 , 2, ... to yield the (noise-free) 
data ..., y(— 1), y(0), y(l), .... Figure 1 shows 
several data samples near the nth knot. Construct a 
computationally efficient interpolating function g(x) 
for which 

1 . g(x) = y(n) when x = n = 0, ±1, ±2, 

2. g(x) is acceptably smooth for all x. 

3. /^[^(x) - g{x)} 2 dx is sufficiently small. 

Condition (1) forces g(x) to interpolate y{x) at its 
samples. If this condition were omitted, we would be 
concerned instead with a version of the standard ap- 
proximation problem. Condition (2) translates into 
a statement about the continuity of g(x) and (per- 
haps) its first and second derivatives. Condition (3) 
is qualitative; that is, what is “sufficiently small”? 
We will address this question quantitatively later. 

There is no loss of generality in assuming that the 
domain of y(x) (and thus #(x)) is all x. We make this 
assumption for mathematical convenience; it enables 
us to ignore the issue of boundary conditions. 

There is a significant loss of generality in making 
the assumption, as is done here, that the sampling 
grid (..., -2, — 1, 0, 1, 2, . . .) is equally spaced. 
Certainly for some functions y(x) the only practical 
way to produce a good interpolating function g(x) is 
to use an unequally spaced sampling grid. However, 
there is an important trade-off involved. As discussed 
in the next section, the assumption of an infinite, 
equally spaced sampling grid enables us to apply 
Fourier (i.e., spectral) analysis techniques to this 
problem and thereby produce valuable insight into 
how to quantify condition (3). 

Many solutions to the interpolation problem 
are available from the numerical analysis, com- 
puter graphics, CAD, and image-processing liter- 
ature. Among the various types of interpolants 
commonly used for data analysis and graphics ap- 
plications are Hermite, cubic, exponential, parabolic, 
rational, B-, beta-, monotonic, and shape-preserving 
splines. For a broad, general introduction to spline 
methods, see reference 1, which contains ready-to-use 
FORTRAN programs. It is important to distinguish 
these spline interpolants from the reconstruction al- 
gorithms used in image processing. The latter have a 
central concept — the interpolation function— around 
which their treatment can be unified. The former, 
however, are more diverse in their derivations and 
mathematical representations, and they are rarely, 
if ever, discussed in terms of their interpolation 
functions. 


The selection of an interpolant is generally a 
highly subjective matter. The choice involves prac- 
tical considerations of computational efficiency and 
subjective evaluations of fairness , that is, the vi- 
sual acceptability of the interpolating function based 
upon its conformance to the data. 

In this paper, condition (3) is used as a criterion 
to measure the accuracy with which a given interpo- 
lating function g reconstructs the original sampled 
function y. The method by which this quantitative 
criterion is used to assess the characteristics of inter- 
polants used for data analysis, CAD, and graphics is 
believed to be new. The criterion is applied to several 
spline interpolants commonly used in data analysis, 
CAD, and graphics as well as to several signal and 
image reconstruction algorithms. 

Although the criterion developed encompasses 
data with arbitrary frequency content, this paper fo- 
cuses primarily on the interpolation of band-limited 
and sufficiently sampled data. However, an exam- 
ple of using the data energy spectrum to select an 
optimal interpolant is presented to indicate how the 
analysis might be applied to tailor a parametric fam- 
ily of interpolants for best performance with a given 
class of data. 


Analysis 

A quantitative analysis of the reconstruction 
properties of interpolants was derived in reference 2 
in an analysis of the effect of sampling and sample- 
scene phase on image reconstruction. The equations 
were then applied in reference 3 to study the image 
reconstruction properties of a particular family of re- 
construction algorithms. In this paper, we show that 
by broadening the interpretation of the basic equa- 
tions in these two references, it is possible to study 
many of the interpolants familiar to the numerical 
analyst and graphics specialist. As a consequence of 
this expanded interpretation, a unified perspective 
is achieved which encompasses both signal and im- 
age reconstruction algorithms and many of the types 
of interpolants used for data analysis, CAD, and 
graphics. 

To begin, we present a summary of the underlying 
equations and concepts from references 2 and 3. The 
main objective of this presentation is to demonstrate 
how to reformulate the mean-square-error criterion 
in a manner which separates the contribution of the 
particular interpolant being used from the contribu- 
tion of the data set being analyzed. Subsequently, 
these equations are applied to a variety of specific 
interpolants. 
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Linear, Shift-Invariant Interpolants 

The analysis of reference 2 was derived for the 
class of interpolants that are both linear and shift in- 
variant. 1 Thus, it is important to clearly understand 
the properties of linear, shift-invariant interpolants. 

An interpolant is linear and shift invariant if and 
only if there exists a function r(x ) such that 

+oo 

0(*)= X y{n)r{x-n) (1) 

n =— oo 

Thus, the output function g(x) may be expressed as 
a linear combination of translated copies of r(x). The 
function r(x) is called the interpolation function for 
the interpolant. It follows from equation (1) that 
each interpolant is uniquely defined by its interpo- 
lation function r(x) and conversely. For a particular 
interpolant, r(x) is the (unique) function which inter- 
polates the data . . . , 0, 0, 1, 0, 0, . . that is, the 
interpolation function r(x) is generated by interpo- 
lating this special data set. The interpolation func- 
tion is fundamental to the quantitative characteri- 
zation of the reconstruction and interpolatory prop- 
erties of its associated interpolant. Reconstruction 
algorithms used for signal and image interpolation 
are usually expressed directly in terms of an inter- 
polation function. For local, linear, shift-invariant 
data and graphics splines, the interpolation function 
can be easily derived in closed form from the defining 
equations by using these data values. 

An interpolant is said to be local if there is an 
integer K such that r(x) = 0 whenever |x| > K. 
The smallest such integer K defines the width of 
the interpolation function, and the interpolant is 
said to be a 2/f-point method. That is, for a local 
interpolant, the infinite sum in equation (1) reduces 
to a finite sum since the terms with \x — n\ > K do 
not contribute; that is, 

[x]+K 

g(z) = X y( n ) r ( x - n ) 

n=[x\— K+l 

where [x] denotes the greatest integer less than or 
equal to x. As an example, the interpolation function 
for linear interpolation is shown in figure 2. Linear 
interpolation is a local method with K = 1. For the 
other local interpolants considered in this paper, K 
is either 2 or 3. 


1 Note, however, that some interpolants which appear in 
the current literature are not linear or shift invariant. An 
example of a nonlinear spline which is outside the scope of 

the analysis in this paper is described subsequently. 


Data, CAD, and graphics splines are frequently 
presented in the literature as a linear combination of 
blending functions with numerical coefficients which 
are functions of the sample values, for example, Her- 
mite interpolation. However, if the numerical coeffi- 
cients are linear combinations of the sample values, 
equation (1) provides an alternate and completely 
equivalent means of carrying out the interpolation. 
The interpretation of equation (1), shown schemati- 
cally in figure 3 for linear interpolation, is that the 
interpolating function, a series of line segments in 
this example, can be built by convolution, that is, as 
a sum of impulse responses , represented by r(x — n), 
of amplitude y(n) centered at the samples x — n. 

Associated with each interpolation function r(x) 
is its reconstruction filter, that is, Fourier transform 
r(is): 

r(v) = f + °° r{x)e~ 2 ^ xu dx (2) 

when j — \f—i. The frequency coordinate v has the 
units of cycles per sample interval and v — 1/2 is 
the Nyquist frequency. Because all the interpolation 
functions r(x) considered in this paper are real and 
even, all the corresponding reconstruction filters r{v) 
are also real and even. 

Mean-Square-Error Criterion 

As discussed previously, we use the mean square 
error, condition (3), 

r+oo 

e 2 = / [y{x) ~ g{x)] 2 dx (3) 

J —OO 

as the criterion which measures how well an interpo- 
lating function g approximates the sampled function 
y. In reference 2 it was shown that if y(x) is band lim- 
ited (i.e., the energy spectrum |y(i/)| 2 is zero for all 
\v\ greater than a cut-off frequency v c ) and if y(x) 
is sufficiently sampled (i.e., v c < 1/2), then equa- 
tion (3) can be transformed into 

/+l/ c 

<?{v)\y{v)\ 2 dv (4) 

- l^c 

In equation (4), 

y{ v) = [ y{x) e ~ 2njx,/ dx (5) 

J — OO 

is the (complex) Fourier transform of the sampled 
function y(x) and 

-hoc 

e 2 (v) = 1 — 2f(i/) + ^ — (6a) 

n=— oo 
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+oo 

e 2 (is) = [1 - r( i/)] 2 + f 2 (i/ - n) (6b) 

n=— oo 

72^0 

Note that, as the notation suggests, e 2 {u) is non- 
negative for all v. Also, if the interpolant is a local 
2/f-point method, it can be shown that an equivalent 
equation for e 2 (y) is 

2K-1 

e 2 (v) = 1 — 2 f{u) + co + 2 ^2 c n cos(27r nv) (7) 

n= 1 

where 

f K 

c n = I r(x)r(n — x)dx (8) 

Jn-K 

If y(x) is band limited but not sufficiently sampled 
(1/2 < u c < oo) or if y(x) is not band limited 
(l/ c = oo), then the mean square error e 2 depends 
upon the location of the function y(x) relative to 
the sampling grid. However, equation (4) remains 
valid in the sense that it represents the average mean 
square error over an ensemble of all possible (random, 
equally likely) shifts of the sampling grid. 

Equation (4) expresses the mean square error as 
an integral of the product of two positive real func- 
tions: e 2 (is) and |y(z/)| 2 . The function e 2 (v) rep- 
resents the contribution of the interpolant, whereas 
the function \y{v)\ 2 represents the contribution of the 
sampled data. Thus, the choice of interpolant has 
been separated from the characteristics of the data. 
Clearly, for a fixed |y(j/)| 2 , reducing e 2 {v) reduces the 
mean square error. In the remainder of this paper, 
we consider the interrelation among r(x), r(z^), and 
e 2 (u) and the potential of minimizing e 2 by reducing 
e 2 (t/). 

Interpolation Methods Investigated 

The interpolants considered in this paper can 
be conveniently grouped into three categories: local 
splines, global splines, and reconstruction filters. The 
spline categories consist of several methods suitable 
for data analysis, CAD, and graphics applications. 
The reconstruction filters were taken from the signal- 
and image-processing literature. Although the inter- 
polants selected were chosen to be representative of 
those which are currently in use or have been studied 
in recent papers, no attempt was made to exhaus- 
tively evaluate all known interpolants. The criterion 
applied here is proposed only as a suitable means 
for estimating the quality of an interpolation with 
linear, shift-invariant interpolants in the context of 
data analysis or signal or image reconstruction. The 


results may not be appropriate for evaluating the per- 
formance of interpolants for many other uses, such 
as rendering or display of data, surfacing, or solid 
modeling. 

Background information on the interpolants along 
with the specifics of their implementation are pre- 
sented in this section. Several of their salient prop- 
erties are summarized for quick reference and easy 
comparison in table I. 

Local Splines 

As discussed previously, the term local spline is 
used to emphasize the fact that the value of the 
interpolating function g at a point x depends only 
upon the values of the 2 K data samples closest to x . 
With the exception of these 2 K samples, the value 
of the interpolant at x is independent of the data. 
For both local splines considered, K = 2, and thus 
only the four nearest sample values contribute to the 
interpolated value. 

Parametric cubic Hermite splines. A three- 
parameter, four-point family of local cubic splines 
with a Hermite basis was introduced in reference 4 
for the purpose of generating “in-between” frames 
from key frames for an automated animation system. 
The three parameters, referred to as “tension,” “con- 
tinuity,” and “bias,” facilitate control of the spline 
shape and enable the user to generate desired ani- 
mation effects. For our purposes, we have chosen to 
set the continuity and bias parameters to zero be- 
cause our experience indicated that nonzero values 
tend to produce distortions unsuitable for data in- 
terpolation. Therefore, only the tension parameter, 
assumed constant for the entire spline (to make the 
spline shift invariant), is used to shape the spline. 
The equations for the resulting one-parameter fam- 
ily of shift-invariant C 1 splines are given in appen- 
dix A. It can be seen from these equations that the 
tension parameter t or equivalently the parameter 
a = — (1 — t) /2, adjusts the values of the spline slopes 
at the end points of the interval containing the point 
x at which an interpolation is desired. From this ob- 
servation, it is clear that although a is unconstrained; 
values of a much outside the range from — 1 to +1 
generally introduce undesirable slopes in g(x) at the 
sample points. When a = — 1 / 2 , one obtains the 
Catmull-Rom spline for which 

y(i + 1) - y{i - 1) y(i)-y{i- 1) , v(i + 1) - V W 

= 2 " 2 + 2 

That is, the slope assigned at the point i is the 
average of the slopes of the left-hand chord y(i) — 
y(i — 1) and right-hand chord y(i + 1) — y(z). 
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Parametric quintic Hermite splines . The one- 
parameter subclass of cubic Hermite splines above is 
identical, except for minor notational changes, to a 
one-parameter family of cubic splines developed inde- 
pendently by Fletcher and McAllister. 2 These splines 
are viewed as data interpolants and are generalized 
in two ways: (1) to two dimensions for generating 
interpolating surfaces and (2) to higher degree. Our 
interest is in this latter extension which gives a two- 
parameter family of fifth degree polynomials. The 
equations for this four-point, local spline are given 
in appendix A. The quintic Hermite basis functions 
H\ (x) to Hq(x) are used; D(i) and D(i + 1) are the 
slopes at the ends of the interval containing the in- 
terpolated point; C(i ) and C(i + 1) are the second 
derivatives. The two parameters t and t! are referred 
to as “tension” parameters by Fletcher and McAllis- 
ter. These parameters control the shape of the inter- 
polant. It can be seen from equations (A5) and (A6) 
that £, or equivalently a = —(1 — t)/ 2, is a weighting 
factor for slope control, whereas t f , or equivalently 
(3 = — (1 — t ! )/2, weights the second derivative terms. 
In this paper a and (3 are assumed constant for the 
entire spline, thus producing a shift-invariant C 2 in- 
terpolant. As with the cubic Hermite splines, values 
of a and (3 out of the range from —1 to +1 introduce 
undesirable slope and curvature terms with conse- 
quent distortion in the interpolating spline. 

Global Splines 

Global splines are fundamentally different from 
local splines in that the value of the interpolating 
function at any point depends upon all the data 
samples. A system of linear algebraic equations 
must be solved in order to construct the interpolant. 
Because of this additional computational effort these 
interpolants are generally used only with relatively 
small data sets. 

Exponential splines . The exponential spline in 
tension (refs. 5 and 6) is a popular form of global 
spline. The tension parameter provides control over 
the spline shape, and the well-known cubic spline is 
obtained as a limiting case of the exponential spline 
with zero tension. Exponential splines are commonly 
used to eliminate the extraneous oscillations which 
are frequently obtained with cubic splines. However, 
exponential splines cannot be evaluated quite as effi- 
ciently as cubic splines, and this has been cited as a 
potential drawback in time-critical applications. 

2 Fletcher, G. Yates; and McAllister, David F.: A Uni- 
fied Approach to Cardinal Spline Surfaces. Unpublished 
research, North Carolina State Univ., under NASA grant 
NAG1-103-S4. 


The construction of an exponential spline requires 
the solution of a tri-diagonal system of linear equa- 
tions which account for the interaction of all the data 
points. The algorithm used in this paper for calcu- 
lating exponential splines is based on the equations 
of reference 6. The tension parameter, designated p 
in reference 6, is held constant along the length of 
the spline in order to preserve shift invariance. 

Nu splines . The nu spline (ref. 7) was devel- 
oped as a piecewise cubic polynomial alternative to 
the exponential spline in tension. This spline is not 
as computationally expensive as exponential splines 
since polynomial functions are involved. Nu splines 
may be generated by using the cubic Hermite basis 
functions in appendix A. The values of the unknown 
slopes at the samples are determined by solving a 
tri-diagonal system of equations which express conti- 
nuity conditions along the spline. For this study, the 
tension parameter, denoted v in reference 7, is con- 
stant along the curve resulting in a C 2 shift-invariant 
interpolant. As in the case of the exponential spline, 
zero tension corresponds to the cubic spline. 

Global interpolation functions . Local splines pos- 
sess spatially limited interpolation functions r(x ) 
whose translation invariance is usually clear from 
the equations defining the spline. The situation is 
quite different for global splines. To generate a global 
spline, one must solve a set of linear equations which 
represent an interaction between all the data values 
and the end conditions. As a result, an interpolated 
value depends upon the end conditions and all the 
sample values, not just upon the values of a few 
neighboring samples. Consequently, global splines 
do not have interpolation functions in the same sense 
as local splines. The function which interpolates the 
data . . . , 0, 0, 1, 0, 0, ... depends upon the location 
of the unit value relative to the end points (i.e., it is 
not shift invariant), the number of samples, and the 
end conditions. Moreover, the interpolation function 
(theoretically) never becomes identically zero. 

Nevertheless, we have verified by numerical exper- 
imentation using from 25 to 100 knots that for both 
exponential and nu splines, interpolation functions 
can be generated which are shift invariant well within 
plotting accuracy (generally, 3-5 significant figures) 
except when centered at a few samples (4-5) near 
the ends. Therefore, as an excellent approximation 
for global interpolants one can define a surrogate in- 
terpolation function which exhibits the essential fea- 
tures of locality and shift invariance and plays a role 
equivalent to the interpolation function for a local 
spline. This function can be calculated in the same 
manner as for a local spline provided that the unit 
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value is not placed too close to an end sample and a 
sufficient number of samples (25 or more) are taken. 
Equation (1) and the subsequent analysis are then 
applicable in a generalized sense, and in this way, we 
are able to include exponential, cubic, and nu splines 
in this study. Because global splines are frequently 
used for curve fitting and data analysis, it is highly 
desirable to establish some basis for comparing their 
reconstruction properties with those of local splines. 

Reconstruction Filters 

Considerable effort has been expended in develop- 
ing efficient and accurate algorithms for interpolat- 
ing and resampling digital signals and images. Since 
the analysis applied here was originally developed in 
an image-processing context, it is natural to include 
several such algorithms in our comparison. These 
algorithms are specified directly in terms of their in- 
terpolation functions, which are given in appendix B. 

Parametric cubic convolution (PCC). PCC is a lo- 
cal, four-point interpolant whose interpolation func- 
tion is given by a one-parameter family of piecewise 
cubic polynomials (eqs. (B2) through (B4)). The 
original analysis of reference 2 was applied to this in- 
terpolant in reference 3 in order to choose an optimal 
value of the parameter a. It was demonstrated that if 
the energy spectrum of the image could be estimated, 
it is possible to choose the parameter to optimize the 
reconstruction properties so as to take this spectrum 
into account. However, for general purpose use, that 
is, in the absence of specific information regarding 
the frequency content of the image, a = — 1/2 was 
found to be the optimum value from the standpoint of 
minimizing the mean square error, provided that the 
data are sufficiently sampled. This choice of a was 
not the value generally referenced in the literature 
and resulted in a reconstruction algorithm generally 
superior to the prevailing standard. 

As indicated in appendix B, the equations for 
the PCC interpolation function are identical to those 
for parametric cubic Hermite splines. Thus, these 
two interpolants are identical. This equivalence does 
not appear to have been observed before and illus- 
trates an advantage to having a unified treatment of 
interpolants. 

Keys 9 cubic. Keys’ cubic is a local, six-point inter- 
polant whose interpolation function is again given by 
piecewise cubic polynomials. However, in contrast to 
PCC, there are no free parameters. The algorithm is 
fourth order convergent (ref. 8) which is the highest 
order which can be achieved with cubic polynomials; 
and therefore, it is expected that this algorithm will 


perform somewhat better than PCC with a = — 1/2. 
However, because this is a six-point interpolant, any 
improvement in performance is achieved at some ad- 
ditional computational expense. 

BAWA cubic . Two 4-point interpolants based on 
piecewise cubic polynomials are presented in refer- 
ence 9. In reference 9, the interpolant whose inter- 
polation function is denoted by ai(£) is identical to 
PCC with a = —1/2. The other interpolant is de- 
fined by a continuous interpolation function with a 
discontinuous first derivative at the sample points 
0, ±1, ±2, a very novel feature. As illustrated in 
reference 9, the associated reconstruction filter has a 
slightly more band-limited spectrum than PCC with 
a = — 1/2. This second interpolant has been included 
in our comparison. 

Quadratic Shape-Preserving Splines 

Because not all local splines are both linear and 
shift invariant, some local interpolants do not fall 
within the scope of the current analysis. For ex- 
ample, parametric cubic and quintic splines lose the 
property of shift invariance if their parameters are 
allowed to vary from one sample interval to the 
next. As the following example illustrates, some in- 
terpolants can also fall outside the scope of the cur- 
rent analysis by failing to be linear. 

The second-degree polynomial “shape-preserving” 
splines described in reference 10 were introduced in 
an effort to preserve the monotonicity and convex- 
ity suggested by a visual inspection of the data. By 
introducing at most one variable knot between sam- 
ples, a quadratic spline is produced which interpo- 
lates the data and preserves shape. This interpolant 
depends upon certain shape characteristics of second- 
degree Bernstein polynomials and the specification 
of appropriate slopes at the knots. The resulting 
spline is C 1 . These splines have also been suggested 
(ref. 11) for shape-preserving interpolation of bivari- 
ate functions on rectangular grids. The examples 
in references 10 and 11 show very smooth interpo- 
lating functions and surfaces free of unwanted wrin- 
kles. There are no variable parameters for controlling 
shape. 

The nonlinearity of this spline arises if the method 
suggested in reference 10 is employed to assign slopes 
at the knots. These slopes M(i) are defined so as 
to ensure that the monotonicity and convexity of the 
data are preserved and are calculated in the following 
manner. At knot £, let S(i) = y(i) — y(i — 1) and 
S{i + 1) = y{i + 1) - y{i) be the slopes of the 
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left-hand and right-hand chords, respectively. Then interpolation function, 


f ° 

M(i) = { 25(i)5(i+1) 
{ 5(0+5(»+l) 


(S(i)S(i+ 1) <0) 
(Otherwise) 


Although the Af(i)’s enter into the definition of the 
spline in a linear manner, they depend upon the 
sample values in a highly nonlinear fashion. The 
spline is local since only four sample points are used 
in an interpolation. However, the resulting spline 
cannot be represented by equation (1), and thus, the 
subsequent analysis does not apply. 


Comparison of Reconstruction Properties of 
Interpolants 

In this section, plots of r(x), r(v), and e(is) are 
presented for the interpolants described previously. 
Examination of the interpolation function r(x) and 
the reconstruction filter, f(i/) provides a qualitative 
indication of the performance of an interpolant. In 
contrast, because of equation (4), the e 2 {u) function 
provides more quantitative information. That is, by 
comparing e(y) functions, judgments can be made 
regarding the minimization of the mean square error 
by the proper choice of interpolant. 

For the data and graphics splines, values of r(x) 
were generated numerically using a 101-point data 
set (100 sample intervals each of unit length) with 
1 for the middle (51st) .data value and all zeros 
otherwise. For the global interpolants, the tri- 
diagonal systems of linear equations given in refer- 
ences 6 and 7 had to be solved. For the local splines, 
the defining equations given in appendix A were used; 
and for the reconstruction filters, the equations in ap- 
pendix B. For each interpolant, the reconstruction fil- 
ter r(i/) was calculated with a fast Fourier transform 
(FFT) routine using 4096 values of r(x) taken over 
an interval that was 100 sample intervals in length. 
Because r(x) is a real, even function of x,r(v) is a 
real, even function of v. Values of the reconstruc- 
tion filter are therefore available for negative v by 
symmetry. Thus, values of r(y) have been calculated 
in steps of 0.01 from v — —20.47 to v — +20.47. 
For \u\ > 20.48, it is assumed that r(u) = 0. From 
equations (6) or (7), it may be seen that e 2 (^) is an 
even function of v. For positive i/, e 2 (u) was cal- 
culated directly from equation (6) with the infinite 
series truncated when \v — n\ > 20.48. 

As an aid in interpreting the results presented 
herein, it is helpful to recall (ref. 12) that if y(x) 
is band limited and sufficiently sampled, then it 
can be reconstructed exactly by using the “ideal” 


sin(7r:r) 
nx 

(M < 0.5) 

(0.5 < M) 

(M < 0.5) 

(0.5 < M) 

as shown in figure 4. 

This reconstruction filter is ideal in the sense that 
if y(x) is band limited and sufficiently sampled, then 
the mean square error e 2 is zero since e 2 {y) — 0 for 
0 < v < 0.5 and |y(^)| 2 = 0 for \v\ > 0.5. However, 
this reconstruction technique cannot be used as the 
basis of a practical interpolation method because of 
the infinite domain of r(z). 

The ideal reconstruction filter defines the desir- 
able characteristics of a practical interpolant to use 
with band limited and sufficiently sampled data: an 
f(i/) which is as flat as possible at unit value for 
v < 0.5 and drops rapidly to 0 for v > 0.5. Or 
equivalently, an e 2 (v) which is as small as possi- 
ble throughout the range 0 < v < 0.5. The lat- 
ter is a quantitative criterion coming directly from 
equation (4). 

Interpolation Function 

The interpolation function r(x) for several of the 
interpolants investigated is shown in figure 5. The 
heavy dots in figure 5 mark the locations of the knots. 
The interpolation functions are all even functions of 
x which take on the special values r(0) = 1 and 

r(±n) = 0 for n = 1, 2, 3, The functions have 

been separated vertically in the figure to facilitate 
comparison. 

A few features deserve mention. The persistent 
oscillation of the ideal interpolant for many sample 
intervals is indicated. The discontinuity in slope of 
the BAWA cubic at the knots is evident. Also, of 
the interpolants studied, the cubic spline has the 
widest interpolation function with an effective K 
of 5; that is, for |x| > 5, |r(x)| < 0.001. As 
the tension is increased from 0 to 50, r(x) for the 
exponential spline changes from that of the cubic 
spline to nearly triangular shape which characterizes 
linear interpolation. 

Figure 6 shows how the shape of r(x) varies 
with a for parametric cubic and quintic Hermite 
interpolation. Recall that a is the slope of r{x) at 
x = 1.0 and note that figure 6(a) also applies to 


for which 


and 
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PCC. Figures 6(a) and 6(b) are quite similar. The 
parameter (3 permits further fine-tuning of the shape 
of r(x) if desired for quintic Hermite interpolation. 

Reconstruction Filter 

Figure 7 presents plots of |r(i/)| on the range 
0 < v < 2.0 for selected interpolants. Keys’ cubic 
(fig. 7(a)), which has the highest order convergence 
that can be obtained with piecewise cubic polyno- 
mials, exemplifies the ideal features of a reconstruc- 
tion filter: flat near v = 0, rapid drop to zero near 
v — 1/2, and only a very weak sideband peaking near 
v = 1.5. 

The reconstruction filter for PCC with a = — 1/2 
(fig. 7(b)) begins to roll off from 1.0 at slightly 
smaller values of v than does Keys’ cubic and has a 
somewhat more gradual decline to 0. There are two 
weak sidebands between 1.0 and 2.0. However, recall 
that PCC uses fewer sample values than Keys’ cubic 
and is therefore a more computationally efficient 
interpolant. 

The reconstruction filter for the BAWA cubic 
(fig. 7(c)) is not as ideal as that for Keys’ cubic. 
However, it is quite similar to PCC with a = — 1/2 
in the range v < 1.0 where it has a somewhat 
narrower spectrum. Note, however, the significant, 
undesirable sideband at u = 1.5. 

Of all the interpolants analyzed, the cubic spline 
displays the best performance for sufficiently sampled 
data as judged by the shape of its reconstruction fil- 
ter. As can be seen from figure 7(d), the reconstruc- 
tion filter for the cubic spline is very flat at v — 0 and 
drops abruptly to zero just beyond v — 1/2. Thus, 
the cubic spline has ideal reconstruction characteris- 
tics superior to both Keys’ cubic and PCC. However, 
because it is a global interpolant, it is considerably 
more computationally intensive than local methods, 
and it is not suitable for reconstructing large amounts 
of data. 

Figures 7(e) and 7(f) show the reconstruction 
filters for two global interpolants. Tightening an 
exponential spline (fig. 7(e)), by adding tension (0 
tension gives the cubic spline) causes the flat section 
at v — 0 to be rounded off and a prominent sideband 
to appear above v = 1.0. These are both undesirable 
features. (As a point of reference, Tension = 50 
for the exponential spline generates nearly linear 
interpolation.) Tensioning a nu spline (fig. 7(f)) 
produces similar effects. 

In summary, of all the interpolants analyzed, the 
one whose reconstruction filter is closest to ideal is 
the cubic spline followed by Keys’ cubic and then the 
BAWA cubic and PCC with a = —1/2. Since the 
cubic spline is a limiting case of the exponential and 
nu splines with 0 tension, there is a small range of 


values of tension just above 0 for which these global 
interpolants also produce more ideal response than 
the local methods. 

The Function e 2 {v) 

The function e 2 {u) appears as part of the in- 
tegrand in the spectral representation of the mean 
square error (eq. (4)). Hence, the shape of e 2 (i/) has a 
direct quantitative influence upon the reconstruction 
performance of an interpolant. Judgments can be 
made regarding the relative mean square error of sev- 
eral interpolants and the interpolant with the least 
mean square error may be identified by comparing 
e 2 {u) functions. Closed form expressions for e 2 (i/) for 
linear interpolation (highly tensioned splines), PCC, 
and the BAWA cubic are given in appendix C. 

Figure 8 shows plots on a linear scale of four 
typical e 2 (v) functions on the range 0 < v < 2.0. 
Recall that e 2 (i /) is an even function of v . The 
curves are very flat near v — 0 but rise rapidly to 
a relative maximum of about 2.0 near v — 1. The 
e 2 (i y) function is almost periodic for large v as can 
be deduced directly from equation (6); that is, as 
\v\ —s ► +oo, \r(i/)\ — ► 0 and the remaining terms 
in equation (6) represent a periodic function with 
period 1. 

Figure 9 is a composite plot of e(y) for several 
interpolants on a semilogarithmic scale for v below 
the Nyquist frequency. This is the range of v which 
contains the significant energy for an adequately 
sampled data set and normally is the only portion 
of the range of integration in equation (4) which 
makes a substantial contribution to the average mean 
square error. Observe from equation (4) that for 
a given data energy spectrum \y(v)\ 2 , the average 
mean square error is minimized by the interpolant 
having the smallest e 2 (^) over this range of v. 

Referring to figure 9, it is again evident that the 
cubic spline is the best choice of all the interpolants 
analyzed. Keys’ cubic is the best of the local meth- 
ods. The BAWA cubic and PCC with a = -1/2 are 
next in line among the local methods. Adding ten- 
sion to either the nu or the exponential spline raises 
the ordinates of the e(^) curve and hence degrades 
the spline’s performance. Small tensions, roughly up 
to about 2.0, result in small deviations from the cubic 
spline. Hence, there is a range of tension for which 
the nu and exponential splines perform better, but at 
greater computational cost, than the local methods. 
This examination of the e[y) function reconfirms the 
conclusions reached by the previous comparison of 
the reconstruction filters. 

Optimized Quintic Hermite Splines 

Parametric quintic Hermite interpolation has two 
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parameters, a and /?. It is of interest to try to select 
these parameters so as to optimize the reconstruction 
performance of this interpolant. In reference 3, 
it was observed that the superior performance of 
PCC with a = —1/2 manifests itself in two ways: 
the reconstruction filter and the e 2 {v) function are 
both very flat near v = 0. Specifically, the Taylor 
series expansions of these functions have the form 
f(z/) = 1.0 + 0(i / 4 ) and e 2 (v) = 0(z^ 6 ) for small v . 
The results of selecting a and (3 so as to duplicate 
this behavior for quintic Hermite splines are briefly 
discussed in this section. 

First, consider choosing the parameters to make 
r(i/) flat at v = 0. The Taylor series expansion of 
r (is) in powers of v was generated from equations (2) 
and (B6) through (B9) using a software package 
for symbolic manipulation. The coefficients of the 
v 2 and i/ 4 terms were then set to zero in order 
to remove the low order terms from the expansion. 
This resulted in two equations for a and (3 with the 
solution a = 117/32 and (3 = 37. The functions r(i /) 
and e(u) for these values of a and /? are shown in 
figure 10. While r(v) is indeed very flat at v — 0, the 
drop to zero at v = 1.0 is delayed to well above the 
Nyquist frequency and a very large sideband becomes 
evident. The e(i /) function is clearly inferior to any 
seen in previous sections. In other words, in this 
case the “design criteria” of choosing the interpolant 
parameters to flatten f (v) near v = 0 was a bad idea. 

Next, an effort was made to choose a and f3 
to make e 2 {v) flat near v — 0. The Taylor series 
expansion of e 2 (v) in powers of v was found up to i/ 6 
and then the coefficients of the v 2 and v 4 terms were 
set to zero and solved for a and (3 giving a = —1/2 
and (3 — — 1. The functions r(i /) and e(v) are shown 
in figure 11. These results produce a spline which is 
certainly quite acceptable for application; that is, the 
design criteria of choosing the interpolant parameters 
to flatten e 2 (i/) near v = 0 proved to be a good 
idea. However, when figure 11(b) and figure 9 are 
carefully plotted on the same scale, it can be seen 
that, while the e(i/) function is close to that for PCC 
with a = —1/2 and that for the BAWA cubic, it 
is actually somewhat above both of these curves. 
Therefore, the optimized quintic Hermite spline is 
not quite as good for minimizing the mean square 
error as either PCC with a — —1/2 or the BAWA 
cubic. 

In conclusion, attempts to optimize the recon- 
struction performance of parametric quintic Hermite 
interpolants by choosing a and (3 to render f(u) and 
e 2 (v) flat near v = 0 have not led to an improved in- 
terpolant. These negative results highlight the need 
to employ more sophisticated optimization criteria in 
order to obtain useful results. 


Optimized Parametric Cubic Interpolants 

In the previous section, an attempt was made to 
optimize the performance of a parametric family of 
interpolants by choosing the parameters to mimic the 
spectral characteristics of the ideal interpolant. In 
this section, we demonstrate that a more comprehen- 
sive optimization of a family of interpolants may be 
performed if \y(v)\ 2 is known. Specifically, we illus- 
trate how PCC can be optimized using an assumed 
model for the data energy spectrum. This example 
illustrates how the theoretical foundation laid earlier 
in the paper may be used to select the parameters in 
a family of interpolants to produce the best recon- 
struction performance. 

Using equation (4) and equation (C7) for e 2 (u) 
for PCC, it can be shown that for data with energy 
spectrum |y(i/)| 2 , the average mean square error of 
the interpolation is 

e 2 = [ eo{u)\y(u )\ 2 du> - 2a ( ei(u)\y{u )\ 2 du 

J — V c J ~V C 

e2M|y(i')| 2 dv (9) 

-v c 


where eo, ei, and are defined in equation (C8) 
and where v c may be oo (if the data are not band 
limited). 

The mean square error is minimum when 
de 2 /da — 0, that is, when 


a ° pt !-uc e 2 {v)\y{v)\ 2 dv 


( 10 ) 


The existence of a minimum is assured since the 
coefficient of a 2 in equation (9) is positive. Thus, 
the mean square error when using parametric cubic 
convolution or parametric cubic Hermite splines is 
least when a — a op t- As shown in reference 3, e\ (y) 
is negative below the Nyquist frequency. Hence, a op t 
is negative for sufficiently sampled data. 

To model the data energy spectrum, we assume 
the model 

l*M' 2 - (II) 

The spectrum, sketched in figure 12, peaks (at a value 
of 2a) at v — 0, then rolls off, and steadily decreases 
to 0 as v — ► oo so as to represent in an ideal way the 
decay of higher spatial frequencies to be expected 
in a typical data energy spectrum. The parameter a 
controls the rate of decay of the higher harmonics and 
the peakedness of the spectrum. For convenience, the 


9 


total energy in the spectrum has been normalized to 
unity: 

f +0C o 

/ \y{v)\ 2 dv=\ 

J —oo 

The unfiltered out-of-band energy (OBE) for any 
data energy spectrum is defined to be the energy 
beyond the Nyquist frequency: 

— 1/2 r+o o 

\y(v)\ 2 dv+ \y(v)\ 2 du 

-oo J 1/2 

For the spectrum in equation (11), 

2 

OBE = 1 tan -1 (7R7) 

7T 

As indicated in figure 12, the OBE measures the 
fraction of energy under that part of the spectrum 
associated with aliasing; large values of OBE are 
associated with significant aliasing. Instead of using 
the explicit parameter a, the implicit parameter OBE 
is used as the independent variable in subsequent 
plots and calculations. The corresponding value of 
a is 


,.l t an[;(l-OBE)]. Trot< ‘ OBE) (12) 

The spectral model in equation (11) exhibits a 
long, slowly decaying tail at high frequencies. In 
practice, much of this tail would have been effec- 
tively removed by appropriate low-pass filtering. We 
allow for that in our spectral model by incorporating 
a function to simulate such low-pass filtering. Specif- 
ically, we assume a three-pole low-pass Butterworth 
filter with its effective cut-off set to the Nyquist fre- 
quency. Our final spectral model, which includes fil- 
tering, is then 


Ml 2 = 


2(7 


1 + (27T(7i/) 2 


l + (2vf\ 


(13) 


The optimal value of a for PCC for this data energy 
spectrum may be calculated from equations (10) and 
(13). In theory, one must set u c = oo to account for 
the aliased part of the data spectrum. In practice, 
u c = 2 is adequate because of the rapid decay of 
equation (13) at high frequencies. Equation (12) is 
used to determine the value of a corresponding to a 
given value of unfiltered OBE. 

The results are shown in figure 13 which presents 
a plot of a op t for out-of-band energies up to 10 per- 
cent (OBE = 0.10). The optimum a is negative over 
most of the range and generally somewhat less than 


the value of —0.5 which is optimum for sufficiently 
sampled (no aliasing) data. A value of —0.7 would be 
an appropriate average to use over the upper 60 per- 
cent of the range of OBE shown. 

It is of interest to compare the absolute mean 
square errors for several of the interpolants studied 
using the spectral model in equation (13). For each, 
the average mean square error has been calculated 
directly from equation (4). The results are shown in 
figure 14 which is a logarithmic plot of e 2 versus un- 
filtered OBE. The curve for linear interpolation also 
represents the results for highly tensioned nu and 
exponential splines. The relative magnitude of the 
mean square errors is seen to be consistent with the 
conclusions drawn earlier in the paper regarding the 
ranking of various interpolants in terms of their re- 
construction performance. That is, the mean square 
error is highest for highly tensioned splines, lowest for 
the cubic spline, with the other interpolants falling 
in between in order as expected. The curve for the 
BAWA cubic is indistinguishable from that for PCC 
with a = —1/2 to the scale used in the plot. Ac- 
tually e 2 for the BAWA is approximately 1 percent 
higher than that for PCC because e 2 (v) for BAWA is 
slightly larger than e 2 (v) for PCC for v > 0.4, indi- 
cating a very slight preference for PCC. It is interest- 
ing to observe that the calculated results in figure 14 
fall on virtually straight parallel lines over the range 
shown implying a relation of the form e 2 ■= A(OBE)^ 
in which A depends on the interpolant but B does 
not. 

Again, recall that PCC, cubic Hermite, and 
BAWA are local four-point interpolants which are 
somewhat more computationally efficient than Keys’ 
six-point cubic and are much more computationally 
efficient than the cubic spline. For the spectral model 
used, these four-point interpolants would be very 
attractive since their mean square errors are only 
slightly larger than those for the more computation- 
ally demanding Keys’ and cubic interpolants. 

Concluding Remarks 

An analysis has been presented which provides 
a quantitative measure of the reconstruction or in- 
terpolation performance of linear, shift- invariant in- 
terpolants. The performance criterion is the mean 
square error of the difference between the sampled 
and reconstructed functions. The analysis is ap- 
plicable to reconstruction algorithms used in signal 
and image processing and to many of the types of 
splines used in numerical analysis, computer-aided 
design (CAD), and computer graphics. When for- 
mulated in the frequency domain, the mean square 
error clearly separates the contribution of the in- 
terpolation method from the contribution of the 
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sampled function and provides a rational basis for 
selecting an “optimal” interpolant; that is, one which 
minimizes the mean square error. Several examples 
of optimized interpolants are presented. 

The analysis has been applied to a collection 
of frequently used reconstruction and interpolation 
techniques: parametric cubic and quintic Hermite 
splines, exponential and nu splines (including the 
special case of the cubic spline), parametric cubic 
convolution, Keys’ fourth-order cubic, and a cubic 
with a discontinuous first derivative (the BAWA cu- 
bic). The performance of these techniques has been 
analyzed in the case in which no a priori knowledge 
of the frequency spectrum of the sampled function is 
assumed and in one specific case where the spectrum 
is assumed known. 

It has been found that with the proper identifica- 
tion of parameters, interpolation with parametric cu- 
bic Hermite splines, developed for computer graphics 
and curve-fitting applications, is identical to recon- 
struction with parametric cubic convolution (PCC), 
a family of reconstruction algorithms used for signal 
and image processing. Of all the interpolants investi- 
gated, it was found that the cubic spline provides the 
best performance. Adding tension to a cubic spline 
in the form of either an exponential or nu spline, gen- 
erally increases the mean square error. Fifth degree 


Hermite splines showed no advantage over third de- 
gree Hermite splines. However, parametric quintic 
splines possess two parameters which might be se- 
lected to tailor an optimal interpolant for a specific 
application. Of the reconstruction algorithms stud- 
ied, Keys’ cubic exhibited the best performance, as 
expected, followed by the very similar BAWA cubic 
with a discontinuous first derivative and cubic Her- 
mite or PCC with a = — 1/2. 

In practical application, these conclusions must 
be weighed against the computational expense. Cu- 
bic, nu, and exponential splines are global inter- 
polants requiring the solution of systems of linear 
equations. Their use with large data sets demands 
considerably more computational resources than the 
local methods — PCC, Keys’ cubic, and BAWA cu- 
bic. Even the choice of Keys’ cubic, which requires 
six samples, in preference to BAWA cubic or PCC, 
which require only four samples, would have to be 
carefully considered for general purpose application 
to high-resolution images because of the additional 
computation. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
March 4, 1987 
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Appendix A 

Equations for Hermite Basis Splines 

This appendix presents the equations for calculating the interpolating function g(x ) for cubic and 
quintic Hermite splines. In the discussion that follows, suppose x lies between i and i + 1 and define 
s — x - i. 

Parametric Cubic Hermite Splines 

The equation for the interpolating function is 

g(x) = y{i)hi(s) + y(i + l)h 2 {s) + D(i)h 3 {s ) + D{i + 1)M S ) (Al) 

where 

D{i) = a[y{i - 1) - y(i + 1)] (A2) 

and a = —(1 — t)/2 where t is the tension parameter used in reference 4. The functions /ij(s) are the 
standard cubic Hermite basis functions defined as 

hi(s) = s 2 (2s - 3) + 1 
h 2 (s ) = s 2 (3 - 2s) 
h^s) = s 2 (s - 2) + s 
h 4 (s) = s 2 (s — 1) 

The h{(s) are defined such that g(i) = y{i), g'(i ) = D(i), g(i + 1) = y(i + 1), and g'(i + 1) = D(i + 1). 
Thus the parameter a weights the slopes of the interpolating function g at the end points i and i + 1, and 
the four neighboring data points y(i — 1), y(i), y(i + 1), and y(i + 2) contribute to the value of g(x). 

Parametric Quintic Hermite Splines 

The equation for the interpolating function is 

g(x) = y(i)H\(s) + y(i + 1 )H 2 (s) + D(i)H 3 (s) + D(i + 1 )H 4 (s) + C(i)H 5 (s) + C(i + 1 )H 6 (s) (A4) 

where 

D{i) = a[y(i - 1) - y(i + 1)] (A5) 

C(i) = -f3[y(i - 1) - 2 y(i) + y{i + 1)] (A6) 

and q = — (1 — t)/2 and (3 = — (1 — t')/2 where t and t' are the two tension parameters. 

The functions H{(s) are the quintic Hermite basis functions: 

H\(s) = 1 — s 3 (6s 2 — 15s + 10) 

H 2 (s) = s 3 (6s 2 - 15s + 10) 

Hs{s) = s — s 3 (3s 2 — 8s + 6) 

H 4 {s) = — s 3 (3s 2 —7s + 4) 

■f^5( 5 ) — — s 2 (s 3 — 3s 2 + 3s — l)/2 
Hq(s) = s 3 (s 2 — 2s + l)/2 

It is easily verified that g(i) = y(i), g'(i ) = D(i), and g"(i) = C(i). Thus, a weights the slopes of g at the 
ends of the interval, and f3 weights the second derivatives. Again, only four data points, y(i - 1), y(i), 
y(i +1), and y(i + 2), enter into the definition of g. 
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Appendix B 

Interpolation Functions 

This appendix presents the equations for the interpolation functions for the local splines and the 
reconstruction algorithms. The functions hi(x) (i = 1, 2, 3, 4) and H{(x) (i = 1, 2, . . . , 6) are the cubic 
and quintic Hermite basis functions defined in appendix A. The equations for the interpolation functions 
are as follows: 

Parametric cubic convolution 


where 


r(x) 


' —ahi{x + 2 ) 
h 2 (x + 1) - ah 3 (x + 1) 
i hi (x) + ahi(x) 
ah 3 {x - 1) 

,0 


(-2 < x < -1) 
(-1 < x < 0) 
(0 < x < +1) 
(+1 < x < +2) 

(2<N) 


= ro(x) + ari(x) 


r 0 (x) = 


(2|x| + l)(|x| — l) 2 (0 < )xj < +1) 

0 (+1 < |x|) 


ri(x) 


|x|uert 2 (|x| — 1) (0 < |x| < 1) 

MN-1)(|x|-2) 2 (1 < |x| < 2) 
,0 (2 < |x|) 


Parametric cubic Hermite splines 

The equations are the same as equations (B2) through (B4). 
Parametric quintic Hermite splines 


r(x) 


' -aH^x + 2) - (3Hq(x + 2) 

H^ix + 1) — aH 3 (x + 1) — f3H$(x + 1) + 2/3Hq(x + 1) 
< Hi{x) + aH A (x) + 2 0H 5 {x) - f3H 6 {x) 
aH 3 (x - 1) - 0H 5 (x - 1) 

,0 


(-2 < x < -1) 
(-1 < x < 0) 
(0 < x < 1) 
(1 < x < 2) 

(2 < N) 


= r 0 (x) + ori(x) + /3r2(x) 


where 


r o(x) = 


— 6|x| 5 + 15|x| 4 

0 


10|x| 3 + 1 (0 < |x| < 1) 

(1 < 1 * 1 ) 


! — 3|x| 5 + 7|x| 4 - 4|x| 3 (0 < |x| < 1) 

-3|x| 5 + 23|x| 4 - 68|x| 3 + 96|x| 2 - 64|x| + 16 (1 < |x| < 2) 
0 (2 < |x|) 


(Bl) 

(B2) 

(B3) 

(B4) 


(B5) 

(B6) 

(B7) 

(B8) 
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r -§|x| 5 + 4|x| 4 - l\xf + |x| 2 (0 < |x| < 1) 

r 2 (x) = < +i|x| 5 - 4|x| 4 + f |x| 3 - 19|x| 2 + 14|x| - 4 (1 < |x| < 2) 

^0 (2 < |x|) 

Keys ? cubic 

' +iM 3 ~ l\ x \ 2 + 1 (0 < |x| < 1) 

= < “i^M 3 + 3|a:( 2 - f§|x| + | (1 < |x| < 2) 

i^M 3 _ |l x l 2 + - i (2 < |x| < 3) 

,0 (3 < |x|) 

BA WA cubic 

' +2;|x| 3 - |x| 2 - ^ |x| + 1 (0 < |x| < 1) 

r(x) = « _1 | x |3 + |a;|2 _ n^i + j (i < | x | < 2) 

,0 (2 < |x|) 


(B9) 


(BIO) 


(Bll) 
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Appendix C 

Some Analytic Expressions for f(i /) and e 2 {u) 

This appendix presents closed form expressions for r(v) and e 2 (v) for several interpolants. These 
formulas were used to validate the numerical procedures used for calculating the values of r(v) and e 2 (is) 
for all interpolants, to confirm the remarkably similar behavior found for PCC and the BAWA cubic below 
the Nyquist frequency, and to check that the results for exponential and nu splines approach those for 
linear interpolation in the limit of large tensions. The expressions for linear interpolation and PCC were 
published previously in references 2 and 3, respectively. The reconstruction filter for the BAWA cubic was 
given in reference 9. The expression for e 2 {u) for the BAWA cubic is new. We use the standard definition 
sine (x) = sin(7r:r)/(7r£). 

Linear interpolation 


r(v) = shur(i/) 
e 2 (v) = 1 — 2 r(v) + 


2 + cos(27 tv) 


(Cl) 

(C2) 


BAWA cubic 


H u ) = 


! (W 


1678 




sinc 4 (i^) 
257 


1680 COS{2niy) ~ TO COS{4niy) + liko C ° S{6ni/) 


(C3) 

(C4) 


PCC and cubic Hermite interpolation 


where 


and 

where 


r(u) = fo(u) + ari(u) 

(C5) 

roM = ^2 [““cfy) sinc(2i/)] 1 

ri(v) = sinc 2 (2i/) — 2 sinc(2i/) — sinc(4i/)j j 

(C6) 

e 2 {v) — eo(^) — 2ae\(v ) + a 2 e 2 (^) 

(C7) 
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eo M = 2- 2 f 0 {u) - — sin 2 (7ri/) 


13 

ei (v) = n(i/) + — sin 2 (2 t tv) 
e 2(v) = j|r sin 2 (27n/) [l + 6sin 2 (7ri/)j 


(C8) 
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Table I. Summary of Properties of Interpolants Analyzed 


Name 

Degree 

Type 

Control 

parameters 

Cubic Hermite 

3 

Local, 4 points 

Slope 

Quintic Hermite 

5 

Local, 4 points 

Slope and curvature 

Exponential 


Global 

Tension 

Cubic spline 

3 

Global 


Nu spline 

3 

Global 

Tension 

PCC 

3 

Local, 4 points 

Slope 

Keys’ cubic 

3 

Local, 6 points 


BAWA cubic 

3 

Local, 4 points 



a 0f both the interpolating and the interpolation functions. 
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Ideal interpolant 




(b) PCC and cubic Hermite, with a — — 1/2, t = 0. 
Figure 7. Reconstruction filters. 
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